clc;
clear;
syms  x y;
z(x,y)=(x^2-2*x)*exp(-x^2-y^2-x*y);
[x0,y0]=meshgrid(-3:.2:3,-2:.2:2);
z0=double(z(x0,y0));
[fx,fy]=gradient(z0);
fx=fx/0.2;
fy=fy/0.2;


figure( 'Name', '引力线图')
contour(x0,y0,z0,30);
hold on;
quiver(x0,y0,-fx,-fy)

zx=diff(z,x);
zx0=double(zx(x0,y0));
zy=diff(z,y);
zy0=double(zy(x0,y0));

disp("-------------------------------------------------")

figure( 'Name', '梯度理论值计算')
subplot(121),
surf(x0,y0,abs(fx-zx0))
axis([-3 3 -2 2 0,0.1])


subplot(122),
surf(x0,y0,abs(fy-zy0))
axis([-3 3 -2 2 0,0.12])
disp("-------------------------------------------------")
figure( 'Name', '误差')
[x1,y1]=meshgrid(-3:.1:3,-2:.1:2);
z1=double(z(x1,y1));
[fx,fy]=gradient(z1);
fx=fx/0.1;
fy=fy/0.1;
z1=double(zx(x1,y1));
z2=double(zy(x1,y1));

subplot(121),surf(x1,y1,abs(fx-z1));
axis([-3 3 -2 2 0,0.1])

subplot(122),surf(x1,y1,abs(fy-z2));
axis([-3 3 -2 2 0,0.1])